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Abstract. We present the results of a three dimensional, locally isothermal, non-self-gravitating SPH code which models 
protoplanetary disks with two fluids: gas and dust. We ran simulations of a 1 M Q star surrounded by a 0.01 M Q disk comprising 
99% gas and 1% dust in mass and extending from 0.5 to ~ 300 AU. The grain size ranges from 10~ 6 m to 10 m for the low 
resolution (~ 25 000 SPH particles) simulations and from 10~ 4 m to 10 cm for the high resolution (~ 160000 SPH particles) 
simulations. Dust grains are slowed down by the sub-Keplerian gas and lose angular momentum, forcing them to migrate 
towards the central star and settle to the midplane. The gas drag efficiency varies according to the grain size, with the larger 
bodies being weakly influenced and following marginally perturbed Keplerian orbits, while smaller grains are strongly coupled 
to the gas. For intermediate sized grains, the drag force decouples the dust and gas, allowing the dust to preferentially migrate 
radially and efficiently settle to the midplane. The resulting dust distributions for each grain size will indicate, when grain 
growth is added, the regions when planets are likely to form. 

Key words, planetary systems: protoplanetary disks - hydrodynamics - methods: numerical 

1 . Introduction migration was fastest for meter sized objects with a velocity of 

~10 4 cms- 1 . 

Protoplanetary disks are by definition where we expect planets | Stenjnskj & Valagea j & ^ added turbulence and 
to form, with micron to millimetre size grains characteristic of 



the interstellar medium collecting and aggregating to form bod- 



grain coagulation and stressed the idea that the dust motion 



,. . ., ,. „ ...... is decoupled from the gas for a given range of particle size. 

les thousands of kilometres in diameter. Grain growth initially „ , , . ... ° . . nA , • 

. . ... . ... t . , They found that grains with sizes 0.1 to 10 cm have inward ra- 

occurs via individual grains sticking via collisions. However, , , r , , ., , 

. ..." . ' , .... , •• dial velocities larger than that of the gas, while larger particles 

once grains reach millimetre size, their collision velocities are , ... 17 : — , ■ „ , 1 J ,»»^4 . , „ , 

ffi . .. . , , . . . , have smaller velocities. Stemnski & Valaaeas ( 1997) modelled 

sufficiently large to shatter the grains upon impact (Jones et al. .... »„ „ , ™ ' , ~ „ . \ • , ,• , 

I. ,[ _ /. , . . , . . ..... a Minimum Mass Solar Nebula (0.24 M Q ) with radial extent 

1996). One way of reducing the relative velocity of colliding , _ . TT . ,. ■ . ■ , 



grains is to increase their number density. Until recently, the 
dust dynamics of protoplanetary disks has mostly been stud- 
ied from the viewpoint of calculating rates of radial migration 

into the central star. Classical work in this field was done by , , . , , 
L T T . . — . .... — |J,„_J, . . , , „ , , settling to the midplane. 
Weidenschilhng ( 1977), who investigated the nature of the drag -j ^ ^ ^ 

force between the dust and gas components of a non-turbulent 



15 AU in which all the solids migrated onto the star, whereas a 
nebula of 0.002 M and extending over 250 AU (which is close 
to our disk parameters) resulted in a distribution of solids close 
to that of our solar system. They did not investigate vertical 



iRice et alJ d2004 studied the concentration of dust in the 



protostellar disk and then derived equations describing the ra- s P iral arms of marginally stable, self-gravitating protoplanetary 

dial motions of individual dust particles. Weidenschilhng con- disks ' ^ followed the evolution of dust test particles added 

eluded that the maximum radial velocity was independent of into their existin 8 SPH code which models § as accretion disks, 

the drag law and was also insensitive to the nebula mass. Radial The test P articles feel the gravitational effects of the star and 

gas disk, as well as gas drag, but do not influence the gas disk. 

Send offprint requests to: L. Barriere-Fouchet They found that the dust density enhancement was significant 
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only for particle sizes between 10 and 100 cm with their neb- 
ula parameters, suggesting that grain growth was possible due 
to the increased density in the dust layer and the decrease in 
the relative velocities of the dust grains to each other. Vertical 
settling was found to be insignificant inside the spiral arms. 

The cas e of vertical settling was investigated by 
Gar aud et alJ (120041) . who analytically derive fluid equations 
for the average motion of particles in a non-turbulent nebula. 
Starting from the momentum equation of a single particle in 
either Stokes or Epstein drag regime (for respectively large and 
small particles), a Boltzmann distribution is used to obtain the 
collective behaviour of dust particles which is consistent as dust 
particles do not interact with each other an d the action of dust 
on the gas is neglected. iGaraud et alJ ll2004 found that small 
particles move smoothly to the midplane, while large bodies 
oscillate about the midplane with decreasing amplitude. 

It should be noted that the grain sizes discussed here are 
larger than those found in the interstellar medium. This is con- 
sistent with the spectral signature recently observed in the disk 
of CQ Tau riTesti et alj EoO.l) suggesting the presence of large 
grains (a few centimetres), as well as evidence of dust process- 
ing and grain growth found in other disks (e.g. iMeeus et alJ 
l2003UAi)aiet all2004 . 

In this article, we present our new three dimensional gas + 
dust code and compare the resu l ts of disk simulations briefly 
to th e work of | We idenschilling| dl977l ) and more extensively 
with lGaraud et alJ d2004 and lGaraud & Linl J2004 . The code 
allows us to follow both the radial migration and vertical set- 
tling of dust in high resolution, and the possibility of including 
additional physical processes such as turbulence, grain evolu- 
tion, and radiative transfer in a simple approximation or de- 
tailed equation of state. It has already been applied in a pre- 
liminary stud y which uses scattered li ght as a diagnosis of dust 
settling (Barriere-Fouchet et al. 2004). 

In a separate project, we incorporate dust physics into a 
parallel N-body + SPH code which calc ulates self-gra vitational 
forces using a tree based data structure (Maddi son et al.l l2003: 
Humble et al .120051) . The data tree adds considerably to the al- 
gorithmic complexity of the code, but means we can consider 
gravitational stability problems. 

In section 2, we present the basic equations for dust hydro- 
dynamics and in section 3 we describe our code. In section 4 
we describe our simulations: the physical hypotheses, units and 
initial state, as well as the results. Finally in section 5 we dis- 
cuss our results in the light of the afore-mentioned references. 



2. Dust dynamics 

A single particle in circular orbit of radius r about a cen- 
tral body of mass M* will move with the Keplerian velocity 
Vk = VGM*/r, where G is the gravitational constant. Gas mov- 
ing in a disk around the same central body will typically orbit 
at a slower velocity due to the radial pressure gradient, which 
solid bodies moving in the gas disk do not experience and so 
orbit at a velocity different than that of the gas. The two phases 
interact via a drag which slows down the dust and forces it to 
migrate inwards to conserve angular momentum. 



Ac cording to Whipple] ill 9721 Il973l) and IWeidenschillind 
d!977l) . small bodies are strongly coupled to the gas and have 
about the same velocity field, whereas large bodies follow 
marginally perturbed Keplerian orbits. Medium size particles 
experience the strongest perturbation to their movement, with 
increased accretion to the central object and vertical settling. 

Before describing the equation of motion, we must explain 
the notations we have adopted for densities. We use the sub- 
scripts d and g to denote dust and gas respectively. We then 
define the intrinsic density (p g , pd) of a fluid and the void frac- 
tion (0 S , 64) as the fraction of the volume occupied by a given 
phase, such that 6d + 8 S — 1 . The density (p e , pd) of a phase in 
a given volume of fluid is then given by 



Pg = °gPg 

Pa = OdPd- 



(1) 



Dust particles have a high intrinsic density (in our case pd = 
1 g cm -3 ), so a given mass of dust occupies a very small volume 
compared to the same mass of gas. Therefore, gas fills almost 
all the volume whereas dust occupies only a small portion of it, 
and: 



Pg ~Pg 
Pd « Pd- 



(2) 



When the only force acting on the dust is drag, the equation 
of motion reads 



dv d Vd - v g mv 
= , with t s = — , 

df ts F D 



(3) 



where Vd is the dust velocity, v g the gas velocity, t the time, f s 
the stopping time and Fd the drag force. The mass of the dust 
grain is 



An , 
m = — p d s 



(4) 



where we have assumed the grains are spherical with radius s. 
We use v = |va - v g | to denote the velocity difference between 
dust and gas phases at a particular point in space. The stop- 
ping time is the time it takes for a dust grain that starts with a 
velocity Vd to reach the gas velocity v g . 

The functional form of the drag force is determined by 
the Reynolds number, Re, and by the ratio of the mean 
free pat h of gas molecules . A, to the radius of the gra ins, s 
(see IWeidenschilhndfl977l [stepinski & Valageaslll996l) . The 
Reynolds number is given by Re = 2sp g v/rj where 77 = p g AC s /2 
is the gas molecular viscosity and C s the sound speed. Thus 
Re = 4sv/AC s . 

If A < As /9, we are in the Stokes regime and the drag is 
due to the wake created by dust particles moving through the 
gas. The expression for the drag force varies according to the 
Reynolds number: 



( 2A Re~ l F for Re < 1 
= \ 2ARe~°- 6 F for 1 < Re < 800 
0.44 F for Re > 800 



with 



F = Ks 2 p g v 2 /2. 



(5) 



(6) 
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If A > As/9 and Re < 1, we are in the Epstein regime and 
the drag is due to thermal agitation and is given by: 



An ~ 
F d = — PgS vC s . 



3. Description of the code 



(7) 



We use the Smoothed Particles Hydrodyn amics (SPH) algo - 

rithm, a Lagrangian technique described by Monaghan ( 1992). 

The SPH equation s and appr oximations have been rigorously 

established bv lBicknelll lll99ll) . 

I I j j 

Our code was based on that of Murray ( 1996), originally 
developed to study tidally unstable accretion disks in cata- 
clysmic variables. We have made several major changes: the 
configuration is for that of a protoplanetary disk rather than 
a binary disk and we have modified the algorithm for finding 
near neighbours so that the code can better hand le variable spa- 
tial res oluti on. Following the work of|Monaghan & Koch arvanl 
i 1995b and lMaddisonl ( 1998), we have added a second "dust" 
phase in a full, self-consistent, fluid approach, contrary to other 
studi es using on ly test particles to describe the dust phase (e.g. 
iRice et all2004l) . We take into account the pressure exerted by 
gas on dust and by dust on gas, as well as the drag force of 
gas on dust. In the results presented in this paper, we did not 
calculate the drag force of dust on gas because it is negligi- 
ble in magnitude and computationally time consuming, but it 
is implemented into the code and can be turned on if required. 
We do not review the SPH meth od as it has b een extensively 
described (see, for example, lMonaghanlll992l). but stress the 
changes made to the code written by Murray ( 1996) in the fol- 
lowing subsections. 



studied bv lOtt & Schnetterl (l2003) and found to be more accu- 
rate. 

In order to calculate the number density and subsequently 
the smoothing length, we need to find all the close neigh- 
bours of a given particle. As our code does not include self 
gravity, rather that usin g a tree code to find neighbours (see 
Hernauist & Katz 1989), we use a less time consuming link 
list. In the original implementation of the code, the cells in 
the linked list were 2h max in size, where h max was the maxi- 
mum value of the smoothing length over the entire simulation 
jMurra\ll 996). With such a cell size, a given particle's neigh- 
bours laid either in the same link cell or in one immediately 
neighbouring it. However, with h varying by as much as a fac- 
tor 1 000, the number of particles in some cells became very 
large and the searches became inefficient. We implemented an 
alternative link-list scheme due to Speith (private communica- 
tion) in which the size of the link cell is chosen to be suffi- 
ciently small so that each cell only contains a few particles. 
This means that, in order to find all the neighbours, the search 
has to be run over a large number of cells. We nevertheless find 
this to be efficient compared to the traditionally used link cell 
algorithm without introducing the programming complexity of 
a tree structure. 

3.3. Drag Force 

For the gas, the equation of motion is given by: 



^ = P ab + M aj + D aj + G„. 



(10) 



P a b is the usual SPH internal pressure term except that each 
fluid is scaled by the volume it occupies, hence 



3.1. Density 

As we are following the evolution of two fluids, we therefore 
have two independent density equations. The gas density is ob- 
tained by summation over only the gas neighbours and the dust 
density is obtained through summation over the dust neigh- 
bours. Using the subscripts a and b to refer to gas particles and 
i and j for dust particles, we find: 



Pi -Xi'"^'n ' 



(8) 



where Wn = W(|r,- - rj\) and W is the cubic spline kernel com- 
monly used in SPH. 

3.2. Smoothing length and link list 

The SPH smoothing length is usually given by 



/ \l/3 
h = h0 [jJ ' 



(9) 



to ensure a roughly constant number of neighbours (in our case 
gas plus dust particles). Since our two types of particles have 
very different masses, to avoid numerical error we therefore 
chose to use the number density n — p/m rather than the mass 
density p in the calculation of h. This approach was extensively 



ab 



z 



m h 



Pcfia PbOb 



Pi 



Pi 



+ H 



nb 



ab, 



(ID 



where Yl a b is the SPH artificial viscosity as described in M urravl 
d!996l) . The volume scaling is ensured by the multiplication by 
the void fraction 9. 

The mixed pressure term, 



t—l n„r>: 



w, 



PaPj 



(12) 



is the pressure exerted on one fluid by the other. 
The drag force term is given by: 



m iK, 



PjPa 



y ja ' ja 



rjaWja, With K a j = ^J^,(13) 



with cr = l/D, D being the number of dimensions, and c a the 
sound speed associated with particle a. The drag exerted by 
dust on gas is negligible and time consuming so we have not 
included it in these computations. 

Finally, G a is the gravity exerted by the central star only, as 
we do not take self gravity into ac count. The derivat ion of the 
M a j and D a j terms can be found in Maddison ( 1998). 

The dust equations are slightly different. The internal pres- 
sure term disappears because the dust fluid is made of large 
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bodies (compared to gas molecules) that hardly ever encounter 
each other. These encounters are better described by an SPH 
shock treatment than by an internal pressure. Thus the dust 
equation of motion is given by: 



dvi 

— - = Ma, + D ib + G h 
df 



where 



V PiPb 
is the mixed pressure term exerted on dust, 



D ib = -<rY J 



m h K ib 
PbPi 



Vib ■ 



ribWib 



(14) 



(15) 



(16) 



is the drag force term exerted on dust, and G, is the gravity 
exerted by the central star. 

We use an implicit scheme which is iterative and hence a 
matrix inversion is not necessary to deriv e the drag force. The 
equations are given bv lMonaghanl Jl997l) and can be found in 
Appendix |A] Note that we use the simpler implicit scheme of 
lMonaghanl ( ll997l) . since our tests of the Tischer algorithm were 
very slow and did not give much better results. 

3.4. Timestep 

A timestep is defined for each physical force: the pressure (<Sf p ), 
gravity {6t g ), and drag (<Sfd) timesteps. The pressure timestep is 
defined so as to verify the Courant co ndition and includes the 
viscosity timestep (see Maddison 1998): 



6t p = min fl 




C a + 0.6^C ab 



(17) 



where subscripts a and b refer to two given SPH particles, c a 
is the sound speed and f a the acceleration (i.e. the net force 
per unit mass) of particle a, and ( is the name we give to the 
SPH artificial viscosity parameter usually called a to avoid any 
confu sion with the Shakura-Sunyaev a ( Sha kura & Sunvaevl 
Il973l) . While 6t v is much larger than 6t g , the gravity subroutine 
is much less time consuming than the pressure subroutine and 
so we therefore use operator splitting. This allows us to enter 
the gravity routine several times while only entering the pres- 
sure routine once. The pressure defines the global timestep 5t, 
such that St = QA6t p , while 6t g gives the number of times the 
gravity routin e is run aft er a single run of the pressure routine, 
following [Murravl lll996h . 

We also use operator splitting for the drag force because 
6t& » 6t g but 5t& can get much smaller than 6t p in case of 
strong drag (which happens close to the central object where 
the gas density is highest and for the smaller grain sizes). We 
choose "by hand" the number of times the drag routine is run 
after a single run of the pressure routine: 5t& decreases with the 
grain size and we find that, down to 1 mm, we just need to run 
it once and then 5 times for 100 fim, 25 times for 10 fim and 
100 times for 1 fim grains. 



3.5. Limitations 

The code treats turbulence in a very limited manner. The dissi- 
pation term in the equation of motion was chosen to reproduce 
the level of dissipation characteristic of the Shakura-Sunyaev 
turbulence model. However, the length scale on which dissipa- 
tion occurs is of the order of the smoothing length and so a full 
turbulent cascade does not develop. 

The resolution of the inner edge of the disk is also lim- 
ited. In the calculations, the disk scale height, which is defined 
by H = C s /£2, where Q = y/GM^/r 3 is the angular veloc- 
ity, has to be larger than the smoothing length h otherwise the 
disk is not resolved. However, H < h at small radii and we 
therefore must set // cot j e = max(H, h) to be able to carry on the 
calculations. It should be noted therefore that the disk is un- 
resolved for r < 6 AU in our high resolution simulations and 
r < 20 - 40 AU in the low resolution simulations. Indeed, our 
disk extends over more than 2 orders of magnitude in radius 
making it difficult to resolve both the small and large radii. The 
solution is not, however, to simply truncate the disk at 6 AU (in 
the high resolution case) because this would result in boundary 
layer problems. Instead we know that we can trust the results 
from 6 AU. Hereafter, what we call the inner parts of the disk 
refer to the regions immediately outside the unresolved parts 
and not the actual disk inner edge. 



4. Simulations 

We model a protoplanetary disk of mass 0.01 M with 1% of 
dust by mass around a 1 M star. We consider dust particles 
ranging from 1 fim to 10 m and run a series of simulations 
that include only one grain size at a time. Our disk is vertically 
isothermal which means that the temperature is constant in the 
vertical direction but follows a radial power law (T cc r~ 3 ^ 4 ). 
We implicitly assume that whatever heating process is acting 
within the disk, the subsequent heat is immediately radiated 
away. 

We do not include the effects of photoevaporation, radiation 
pressure, Poynting Robertson flux, magnetic fields or grain co- 
agulation, sublimation or condensation. 

We scale our model so as to have numbers close to unity. 
The mass unit is one solar mass, the length unit 100 AU and the 
time unit is the orbital period of a particle at 100 AU around a 
1 M object divided by 2k (1000 yr/27r). This choice of units 
leads to a gravitational constant set to 1 . 

We choose an initial state in near equilibrium conditions. In 
the following simulations, we let a gaseous non dusty disk relax 
from an initial distribution given by X = Zor -3 ^ 2 and H — C s /O. 
The initial velocity is Keplerian, given by vj = y/GM+fr. For 
the power laws of the temperature and initial surface density, 
we take the pa rameters of the Minimum Mass Solar Nebula 
llHavashietalll985l) . 

Starting from that initial distribution, we let the disk relax 
for ~ 8 000 years, i.e. 8 orbits at 100 AU, which allows the 
pressure and artificial viscosity time to smooth out the velocity 
field. The gas velocity becomes slightly sub-Keplerian because 
of the pressure gradient. Once the gas disk has relaxed, we then 
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X (AUJ 



Fig. 1. Density contours for dust particles for the disk seen face on at the end of the low resolution simulations, each quadrant 
representing a different size of particles. The vertical bar gives logp in g crrT 3 . 
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add the dust particles on top of the gas particles with the same We checked that the Reynolds number is always less than 
velocity and let the disk evolve for a further 8 000 years. unity since the velocity difference between dust and gas parti- 

cles stays subsonic. The mean free path is given by 
1 m H , 

A= *— (18) 
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Table 1. Simulation parameters 





= 1.0 M 


.V 


= 10,1,10- 1 ,10- 2 ,10- 3 ,10- 4 ,10- 5 ,10-* m 


M m 


= 0.01 Mj 


c s 


= C r- 3/s , C = 0.1 code units 


A/dust 


= 0.01 M gas 


2 




-^disk 


= 300 AU 


«SPH 


= f = 0.1 


Pi 


= 1.0 gcirr 3 




= 0.0 



The right hand sid e is obtained u sing the approximations made 
by ISteoinski & Valaeeasl Jl99d) . that hydrogen is considered 
as molecular and the cross section is assumed to be identical 
to the geometrical section. A decreases with increasing density, 
so if the criterion for the Epstein regime, A > 4s/9, is met 
in the densest regions (i.e. the centre of the disk) then it will 
be met everywhere. Because our disk is not massive (Md; s k = 
0.0 1M Q ) and is radially extended (300 AU), the highest gas 
volumetric density we get is of about 5 10~ n gem -3 , giving 
9 A/4 > 11.5 m. Therefore, for 10~ 6 m < s < 10 m, we are in 
the Epstein regime, and 



h = 



and the equation of motion for the dust is given by 

dva C s p g 
df 



-T- = (Vd -Vg). 



spd 



(19) 



(20) 



We studied a large range of dust grains sizes at low resolu- 
tion (~ 25 000 SPH particles) to see the effect of grain size on 
the dust settling morphology. The simulations proved to have 
evolved for long enough to see the most striking features of 
the dust settling. We then focused on intermediate size grains 
(10 cm to 100 //m) where the drag is most efficient and ran high 
resolution simulations (~ 160000 SPH particles). The results 
are shown on Figs. Qand[2| for the low resolution, and Tabled 
lists the simulation parameters. 

A high resolution single grain size computation requires an 
average of 15000 CPU hours shared over 64 processors with 
OpenMP parallelisation on the SGI Origin 3 800 of the French 
national computing centre CINES based in Montpellier. As the 
grain size decreases, the drag timestep gets smaller and the 
computation time increases. The 100 /mi simulation required 
22400 CPU hours. 

5. Discussion 

Our simulations show that the dust in th e solar n ebula be haves 
qualitatively in the way of the analysis o fWhiDpie1 lfl972ll 19731) 
and IWeidenschillind Jl977h . For example, the larger bodies 
(s > 1 m) are weakly coupled to the gas and follow marginally 
perturbed Keplerian orbits and the dust disk is flared except 
near the centre where the gas drag is the most efficient because 
of the high gas density. On the other hand, the smaller bodies 
(s < 10 yum) are so strongly coupled to the gas that they follow 
its motion and again the dust disk is flared. 

But for intermediate size particles (100 /mi < s < 10 cm), 
the drag strongly influences the dust dynamics and induces a 
motion that is very different from that of the gas. We first see a 
rapid stopping phase, where the dust settles into a region where 




200 



Fig. 3. /< = 1 contours (see text) at the end of the high- 
resolution simulations for 100 //m, 1 mm and 1 cm grains. For 
the 10 cm particles, this value was not reached and the contour 
is drawn for^u = 0.5 instead. 



the drag force dominates. This region is shown in Fig.[3]as the 
interior o f the u = 1 contour s. This parameter fi was charac- 
terised bv lGaraud et alJ (120041) as the ratio of the orbital time to 
the stopping time: 

J_PgQ 



(21) 



When fx » 1, the Epstein drag is very strong. Once the larger 
grains reach the /i > 1 region, they experience rapid settling 
to the midplane as well as strong radial accretion and fall in 
the central star. The dust layer therefore gets very thin at small 
radii. For the smaller grains, the settling to the midplane is so 
efficient that radial migration cannot respond fast enough and 
there is a particle pileup and the dust layer gets thicker again. 
Unfortunately, we cannot quantitatively predict the thickness 
of the dust layer with these simulations, as dispersion in the 
velocities of individual SPH particles limits our resolution. 

Figure^also shows the transition between regimes of weak 
and strong vertical settling, as well as weak and strong radial 
migration. The 10 m particles do not show any particular in- 
crease in density because the settling is weak, and their distri- 
bution stays very close to their initial configuration, which is 
that of the gas (as one can see in the edge-on plots in Fig.|2j. 
Then for 1 m particles, the settling becomes efficient in the in- 
ner regions of the disk, but the radial migration is still weak. 
The particles therefore pile up towards the inner edge of the 
disk and the density increases. For 10 cm particles, the settling 
is even stronger and extends to larger radii, but then the radial 
migration becomes efficient and the particles are lost to the star 
and hence the density enhancement is weaker than for larger 
grains. For 1 cm particles, the settling is now so efficient that 
it extends over even larger radii and despite efficient accretion, 
too many particles arrive at the inner part of the disk at the same 
time and thus they pile up and we see a strong density enhance- 
ment. For 1 mm particles, the combined effects of settling and 
migration have proven so efficient that the outer parts of the 
disk are depleted of dust. Then the efficiency of the vertical 
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10 100 
r (AU) 



Fig. 4. Surface density £ (top panel) and volume density p (bot- 
tom panel) as a function of r at the end of the high-resolution 
simulations. In each case, the dust density is multiplied by 100 
to better view the dust density increase. 

settling decreases together with the radial migration, up to the 
point where particles are completely coupled to the gas as for 
1 fim particles, for which the distribution is undistinguishable 
from their inital configuration. 

In Fig. |2 one can see a large increase in the volume den- 
sity ratio between dust and gas, rising from 10~ 2 to 10" 1 for 
some radii. This increase is not as striking when we look at 
the surface density because the effect of the settling is hidden 
by the summation over the vertical scale height of the disk. 
Therefore 2D simulations will underestimate the density en- 
hancement. The meaningful quantity to study is in fact the vol- 
umetric density, which has the potential to sufficiently damp 
the dust velocities to allow the grains to stick together via col- 
lisions instead of being shattered. It should be noted, however, 
that while we do not consider the drag of the dust on the gas 
because it is negligible when the dust to gas ratio is small, as 
this ratio approach one the dust could drag the gas with it into 
high density regions, thus resulting in a self regulation of the 
dust to gas density ratio. 

iHumble eTaiT d2005h show that, for values typical of the 
solar nebula, the distribution of solid particles evolves to an 
essentially stationary state and that the final radius of the dust 
disk can b e characterised in terms of the particle size. 

As in IWeidenschillind Jl977l) . we find that the radial mi- 
gration rate is highest for a given grain size and drops quickly 
for larger and smaller particle sizes. We fin d this maximum 
is reac hed for 1 cm size particles, whereas Weidenschilling 
il 19771) found it to be 1 m. This discrepancy comes from 
the difference in the p arameters we use for the nebula and 
IWeid enschillingl i 19771) showed that the gas density and dust 
intrinsic density have an impact on this maximum grain size. 



10~ 16 ii i i i | 




z (AU) 

Fig. 5. Volume density versus z at r = 150 AU for the 100 fim 
dust at different times from dust injection (lightest curve) until 
the end of the high-resolution simulations (darkest curve). 

To check the vertical evolution, we first compare our results 
with the self-similar behaviour of the dus t phase during set- 
tling as described bv lGaraud & Linl J2004) and shown in their 
figures lb and 2. Despite the noise due to the SPH technique, 
our results indeed show the same kind of behaviour, shown in 
Fig. [5] Note that our density increase is orders of magnitudes 
smaller than theirs because the computation is too time con- 
suming to be run fo r as long an evolutionary time as given in 
lGaraud&Lin1ll2004l) . 

We next validate the formula derived by iGaraud et alJ 
J2004) in the Epstein regime and reproduce their figure 5 
with our conditions (see Fig. [6ji. Starting from the trajectory 
of a single particle and then using a Boltzmann averaging, 
IGaraud et alJ J2004I) derived an analytical formula for the ver- 
tical velocity 

V z = --, (22) 

with // defined in Eq. Mil . To calculate the value of fi and then 
of v z , we use the gas density that comes out of the simulations 
and the sound speed prescription that we use in our code. We 
then compare this computed value of v z with the one taken di- 
rectly from the simulations in Fig. [6] 

The proportionality between vertical velocity and vertical 
position is most striking at the very beginning of the simula- 
tion when the vertical extension is maximal and the particle ve- 
locities towards the midplane are largest. The plots are shown 
400 yrs (0.4 orbits at 100 AU) after the addition of dust. 

We can easily see the vertical stratification of the velocity 
for the 1 mm grains. It is less obvious for the 100 fj.m grains 
because the dust is more strongly coupled to the gas and set- 
tles more slowly towards the midplane. Furthermore, the noise 
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from the simulation makes the structure more difficult to de- 
tect. At the beginning of the simulations for the 10 cm and 1 cm 
particles, the grains oscillate about the midplane following per- 
turbed Keplerian orbits before the gas da mps these os cillations 
efficiently as shown on figure 2 of Gara ud et alJ (12004). 

Finally, looking at the sequence of edge-on images in 
Fig. El we see q ualitatively the behaviour discussed by 
iGaraud et alJ J2004I) that at a given height the region will be 
depleted of intermediate size particles. 

6. Conclusion 

6.1. Summary 

In this article, we have shown the importance of the vertical 
settling of intermediate size grains (10 cm to 100 yum) in a 
protoplanetary disk with a central object of M* = 1.0 M Q , a 
gas disk mass of M gas = 0.01 M* and a dust disk mass of 
Mdust = 0.01 M g 

as- 

We find that the dust distribution is completely different 
from the gas and the general approximation that the dust to 
gas ratio is constant throughout the disk is wrong for grains in 
this size range. This has an impact on the interpretation of ob- 
servations made in the millimetre domain, for instance, as the 
opacity of millimetre size grains dominates. It is thus important 
to understand the specific behaviour of the dust to avoid mis- 
interpretations such as an underestimate or overestimate of the 
gas density. 

Our simul ation results compare w ell with the radial mi- 
gration seen in Weidenschillin 2 i 19771) . the analyti cal work of 
IGaraud et alJ J2004I) and the vertical settling of IGaraud & LinI 



6.2. Perspectives 

The lShakura & Sunvaevl i ll 9731) turbulence modelling has been 
extensively used in the domain of protoplanetary disks and 
has proven to give consistent results. It should, however, be 
used cautiously and ideally a better prescription for turbulence 
should be used. Indeed, turbulence will likely reduce the dust 
settling velocities and lessen the densi ty enhancements we ob- 
serve in our simulations (see, e.g., IWeidenschilling & Cuzzil 
119931) . The magneto-rotational in stabili ty has been consistently 
described by Balbus & Hawlevl dl991 ) and global numerical 
simulatio ns for this instability are becoming more accurate 
jFromanget al.l l2004). It is now possible to have a consistent 
description of this kind of turbulence in a numeri cal simula- 
tion an d a first step towards this goal was given by Monaghan 
J2002I) . who derived an SPH prescription of turbulence. We 
plan to incorpate SPH turbulence into our code in order to in- 
vestigate how turbulence affects the dust settling. 

While the locally isothermal approximation is roughly 
consistent with the temperature distribution in a proto- 
planetary disk, a lot of effort is currently going into the 
domain of radiative transfe r in protoplanetary disks (e.g. 
iDullemond & Dominikll2004T) . suggesting that the energetics 
in such an object deserve a better treatment. Unfortunately, 
including full radiative transfer in our SPH code is not fea- 



sible because of the increased computation it would require. 
Nevertheless, we will be able to achieve a better description by 
implementing an adiabatic equation of state with cooling func- 
tions as one of the improvements we plan to add in our code. 

This consistent description of the dust distribution in a pro- 
toplanetary disk is a first step towards planet formation. By 
implementing a consistent treatment of the grain growth, co- 
agulation, and shattering, we will be able to better understand 
planetesimal formation and their distribution in the disk. 
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Appendix A: Implicit schemes for the drag force 

All the following equations are derived from Monaghan ( 1997). 
The momentum equations for gas and dust in Epstein regime 
read 



dv g K 

-IT = ~ — (Vg-Vd) 

df p g 
dv d K 

-7- = ~ — (Vd-Vg), 

df p d 



with 



K 



CsPgPd 

pas 



(A.l) 
(A.2) 

(A.3) 



and become in SPH notations (indices a and b refer to gas 
particles and i and j to dust particles and for a vector q, 
q a j = q a ~ Qj)- 



= (T > TO; 



PjPa 



dvj v-^ K aj 
— = cr > m fl — 
dt ^ p a pj 



. r 2 . + n 2 I 

V aj l J 



r ajW a j 



r jaWja- 



(A.4) 



(A.5) 



Then, it is possible to turn it to an implicit scheme without 
having to invert a matrix. We have 



dv fl v 



dv 



— = 2j m aSaj(y a j ■ r aJ )r ja , 



dt 
with 



K a j 
r 

PaPj 



( W aJ " 

r 2 . + n 2 

aj I ! 



(A.7) 



(A.8) 



With a timestep St and v° and v\ respectively the old and new 
values of the velocity: 



j 



(A.9) 
(A. 10) 



Now, we will just consider one dust-gas particle pair at a time: 
'«.,,"«,/ (A.ll) 

(A.12) 



v x a =vl-StmjS a j{vl r r a j)r, 



v)=v°.-6tm a s aJ (v 1 aj -r a j)r Ja . 



Then taking the difference of the previous two equations and 
taking the scalar product with r a j, we get 



V aj r "J 



1 +5t(m a +mj)s aj r 2 a ' 



and 



<0 8tmjs aj r a j(v aj ■ r aj ) 
1 +6t(m a +mj)s aj r 2 a 



v 1 = v° - 

j J 



6tm a s aj r ja (v aj ■ r aj ) 
1 + 6t(m a + rrij)s a jr 



(A.13) 



(A. 14) 



(A.15) 



j" 



The drag of dust on gas is very small, and to save computing 
time, we neglect it. So, at the end, only dust experiences drag 
with the following expression: 

ZStm a s a ;r ijy . ■ r„i) 
(A.16) 
. 1 + 8t(m a + mj)s aj r 2 ja 



A.1. Tischer implicit scheme 

Now, we subdivide each timestep into two, and we note v the 
velocity at A = St/2: 



v a = v° - Ar aj mjS a j(0.6v aj ■ r aj + 0.4v^ ■ r aj ), 



(A. 17) 



vj =v°j- Ar aJ mjS a j(Q.6v aJ ■ r aj + 0.4v° . ■ r aj ), (A. 1 8) 

and 



v, 1 , = 1.4v fl - 0.4v° - 0.6Ar a jmjSaj(ylj ■ r aj ), 
v) = lAvj - 0Av°j - 0.6Ar aj m a s a j(vlj ■ r aj ). 



(A. 19) 



(A.20) 



(A.6) For simpler notations, we define A = Ar 2 j(nij + m a )s a j and 



B = Ar a jin a s a j. Then 



_ 1 - 0.4A 
Vaj ' r ° j ~ 1+0.6A V ^' ' r " h 



1 - 0.8A o 



v >J' r 'l- ( l + .6A) 2V ^"" fly ' 



and 



v' = v" + Br, 



2 + 0.36A 



^"' aj (l+0.6A) 2V "j' raj - 



(A.21) 
(A.22) 

(A.23) 



